/*-------------------------------------------------------------------------*/
/**
  @file		interface_metis.c
  @author	Renaud Wanschoor
  @date	july 2004
  @version	$Revision: 2.4 $
  @brief	interface between metis and plib functions

  A few functions are used in common between different partitionners, those in
  md_io.c are a example, also pagridplot is able to use data from PAGRID, JOSTLE
  METIS and it uses interface to be able to use tool specific data
  
*/
/*--------------------------------------------------------------------------*/

#include "plib.h"

#define RandomInRange(u) ((int)(drand48()*((double)(u))))
#define RandomInRangeFast(u) ((rand()>>3)%(u))

#define SWAP(a, b, tmp)  \
                 do {(tmp) = (a); (a) = (b); (b) = (tmp);} while(0)

#ifdef METIS
static GraphType* __Init(void);
static void __Allocate(GraphType* g);
static GraphType* __Read(FILE *fp);
static void __Free(GraphType* graph);
static GraphType* __Buildcoarse(GraphType* graph);

GraphType* __Init(void) 
{
  GraphType* graph;
  graph = U_malloc(1,GraphType,"Metis_Init error code 1\n");
  
  graph->gdata = graph->rdata = NULL;

  graph->nvtxs = graph->nedges = -1;

  graph->xadj = graph->vwgt = graph->adjncy = graph->adjwgt = NULL;
  graph->adjwgtsum = NULL;
  graph->cmap = NULL;

  graph->ncon = -1;
  graph->nvwgt = NULL;

  graph->vsize = NULL;
  graph->where = NULL;

  graph->npwgts = NULL;

  graph->coarser = graph->finer = NULL;

  return graph;
}

void __Allocate(GraphType* g){
  int *data;

  data = U_malloc(5*g->nvtxs + 2*g->nedges + 1,idxtype, "Metis_Allocate error code 2\n");

  g->gdata     = data;
  g->xadj      = data;
  g->vwgt      = data + (g->nvtxs + 1);
  g->vsize     = data + (2*g->nvtxs + 1);
  g->adjncy    = data + (3*g->nvtxs + 1);
  g->adjwgt    = data + (3*g->nvtxs + g->nedges + 1);
  g->adjwgtsum = data + (3*g->nvtxs + 2*g->nedges + 1);
  g->cmap      = data + (4*g->nvtxs + 2*g->nedges + 1);

}

GraphType* __Read(FILE *fpin){
  int i, j, k, l, fmt, readew, readvw, ncon, edge, ewgt=1, wgtflag;
  idxtype *xadj, *adjncy, *vwgt, *adjwgt;
  char *line, *oldstr, *newstr;
  GraphType* graph;

  graph = __Init();

  line = U_malloc(MAXLINE+1,char,"MD_MetisReadGraph error code 1\n");

  do {
    fgets(line, MAXLINE, fpin);
  } while (line[0] == '%' && !feof(fpin));

  if (feof(fpin)) {
    graph->nvtxs = 0;
    U_free("MD_MetisReadGraph error code 2\n",&line,LEND);
    return graph;
  }

  fmt = ncon = 0;
  sscanf(line, "%d %d %d %d", &(graph->nvtxs), &(graph->nedges), &fmt, &ncon);

  readew = (fmt%10 > 0);
  readvw = ((fmt/10)%10 > 0);
  if (fmt >= 100) {
    printf("Cannot read this type of file format!");
    exit(0);
  }


  wgtflag = 0;
  if (readew)
    wgtflag += 1;
  if (readvw)
    wgtflag += 2;

  if (ncon > 0 && !readvw) {
    printf("------------------------------------------------------------------------------\n");
    printf("***  I detected an error in your input file  ***\n\n");
    printf("You specified ncon=%d, but the fmt parameter does not specify vertex weights\n", ncon);
    printf("Make sure that the fmt parameter is set to either 10 or 11.\n");
    printf("------------------------------------------------------------------------------\n");
    exit(0);
  }

  graph->nedges *=2;
  ncon = graph->ncon = (ncon == 0 ? 1 : ncon);

  /*printf("%d %d %d %d %d [%d %d]\n", fmt, fmt%10, (fmt/10)%10, ncon, graph->ncon, readew, readvw);*/

  if (graph->nvtxs > MAXIDX){
    fprintf(stderr,"\nThe matrix is too big: %d [%d %d]\n", graph->nvtxs, MAXIDX, sizeof(idxtype));
    return NULL;
  }

  xadj = graph->xadj = U_malloc(graph->nvtxs+1,idxtype, "MD_MetisReadGraph error code 3\n");
  for(j=0;j<graph->nvtxs+1;j++)
    xadj[j] = 0;
  adjncy = graph->adjncy = U_malloc(graph->nedges,idxtype,"MD_MetisReadGraph error code 4\n");

  vwgt = graph->vwgt = (readvw ? U_malloc(ncon*graph->nvtxs,idxtype,"MD_MetisReadGraph error code 5\n") : NULL);
  adjwgt = graph->adjwgt = (readew ? U_malloc(graph->nedges,idxtype,"MD_MetisReadGraph error code 6\n") : NULL);
  graph->cmap = U_malloc(graph->nvtxs,idxtype,"MD_MetisReadGraph error code 7\n");

  /* Start reading the graph file */
  for (xadj[0]=0, k=0, i=0; i<graph->nvtxs; i++) {
    do {
      fgets(line, MAXLINE, fpin);
    } while (line[0] == '%' && !feof(fpin));
    oldstr = line;
    newstr = NULL;

    if (strlen(line) == MAXLINE){
      fprintf(stderr,"\nBuffer for fgets not big enough!\n");
      return NULL;
    }

    if (readvw) {
      for (l=0; l<ncon; l++) {
        vwgt[i*ncon+l] = (int)strtol(oldstr, &newstr, 10);
        oldstr = newstr;
      }
    }

    for (;;) {
      edge = (int)strtol(oldstr, &newstr, 10) -1;
      oldstr = newstr;

      if (readew) {
        ewgt = (int)strtol(oldstr, &newstr, 10);
        oldstr = newstr;
      }

      if (edge < 0)
        break;

      adjncy[k] = edge;
      if (readew) 
        adjwgt[k] = ewgt;
      k++;
    } 
    xadj[i+1] = k;
  }

  if (k != graph->nedges) {
    printf("------------------------------------------------------------------------------\n");
    printf("***  I detected an error in your input file  ***\n\n");
    printf("In the first line of the file, you specified that the graph contained\n%d edges. However, I only found %d edges in the file.\n", graph->nedges/2, k/2);
    if (2*k == graph->nedges) {
      printf("\n *> I detected that you specified twice the number of edges that you have in\n");
      printf("    the file. Remember that the number of edges specified in the first line\n");
      printf("    counts each edge between vertices v and u only once.\n\n");
    }
    printf("Please specify the correct number of edges in the first line of the file.\n");
    printf("------------------------------------------------------------------------------\n");
    exit(0);
  }

  U_free("MD_MetisReadGraph error code 3\n",&line,LEND);
  return graph;
}

void __Free(GraphType* graph){
  if(graph == NULL)
    return;

  if(graph->gdata != NULL)
    U_free("MD_MetisFreeGraph error code 1\n",&(graph->gdata),LEND);
  else{
    U_free("MD_MetisFreeGraph error code 2\n",&(graph->xadj),&(graph->vwgt),
	    &(graph->vsize),&(graph->adjncy),&(graph->adjwgt),
	    &(graph->adjwgtsum),&(graph->cmap),&(graph->nvwgt),LEND);
  }
  if(graph->rdata != NULL){
    U_free("MD_MetisFreeGraph error code 3\n",&(graph->rdata),LEND);
    graph->where = NULL;
  }
  if(graph->where != NULL){
    U_free("MD_MetisFreeGraph error code 4\n",&(graph->where),LEND);
  }
  U_free("MD_MetisFreeGraph error code 5\n",&graph,LEND);
}


/**
 * \brief translate a specific partitionner application graph into a ::PlotGraphType graph
 *
 * The returned graph must be freed using PLIBG_Free()
 * \param g application graph
 * \return ::PlotGraphType graph
 */
PlotGraphType* Metis_TransGraph(void* g){
  int i;
  GraphType* graph = g;
  PlotGraphType* pgraph;

  pgraph = PLIBG_Init();

  pgraph->nvtxs       = graph->nvtxs;
  pgraph->nedges      = graph->nedges;

  pgraph->adjidx      = graph->xadj;
  pgraph->adjncy      = graph->adjncy;
  if(graph->vwgt){
    pgraph->vwgt      = graph->vwgt;
    pgraph->adjwgt    = graph->adjwgt;
    pgraph->adjwgtsum = graph->adjwgtsum;
  }else{
    pgraph->vwgt      = U_mallocset(graph->nvtxs,int,1,"Metis_TransGraph error code 1\n");
    pgraph->adjwgt    = U_mallocset(graph->nedges,int,1,"Metis_TransGraph error code 2\n");
    pgraph->adjwgtsum = U_malloc(graph->nvtxs,int,"Metis_TransGraph error code 3\n");
    for(i=0;i<graph->nvtxs;i++){
      pgraph->adjwgtsum[i] = graph->xadj[i+1] - graph->xadj[i];
    }
  }
  pgraph->map         = graph->cmap;
  pgraph->pro         = graph->where;

  return pgraph;
}

/**
 * \brief translate a specific partitionner grid into a ::PlotGridType grid
 *
 * The return grid must be freed using PLIBN_Free()
 * \param g computational grid
 * \return ::PlotGridType grid
 */
PlotGridType* Metis_TransGrid(void* g){
  GridType* grid = g;
  PlotGridType* pgrid;

  pgrid = PLIBN_Init();

  pgrid->npro        = grid->npro;
  pgrid->nedges      = grid->nedges;
  pgrid->ratio       = grid->ratio;
  pgrid->exec        = grid->exec;
  
  pgrid->r_from      = grid->r_from;
  pgrid->r_to        = grid->r_to;
  pgrid->r_val       = grid->r_val;

  pgrid->netcost     = grid->realnetcost;
  pgrid->proratio    = grid->proratio;
  pgrid->proportion  = grid->proportion;
  pgrid->prospeed    = grid->prospeed;
  pgrid->penalty     = grid->penalty;

  return pgrid;
}
/**
 * \brief Serialize a graph into a file
 *
 * \param graph application graph to be serialized
 * \param filename file in which the graph is to be serialized
 * \param outresult parameter for multiple out put (see IO_Args())
 * \param id id of the multiple output file (see IO_Args())
 * \param F_graph string that points to the graph file
 * \param F_grid string that points to the grid file
 * \param CType IoCtrl::CType
 * \param WType IoCtrl::WType
 * \param RType IoCtrl::RType
 * \param CSeed IoCtrl::CSeed
 * \param ISeed IoCtrl::ISeed
 * \param RSeed IoCtrl::RSeed
 */
void Metis_Serialize(GraphType *graph,char *filename,
		      int outresult, int id,
		      char* F_graph, char* F_grid,
		      int CType, int WType, int RType,
		      int CSeed, int ISeed, int RSeed){
  int i;
  int cnvtxs, nvtxs, *match, *map, count = 0;
  GraphType* cgraph;
  FILE* fp;

  while(graph->coarser != NULL)
    graph = graph->coarser;

  cgraph = graph;

  while(graph->finer != NULL){
    count++;
    graph = graph->finer;
  }

  fp = PLIB_SerializeInit(TOOL_METIS,count,filename,outresult,id,F_graph,F_grid,CType,WType,RType,CSeed,ISeed,RSeed);

  if(count>0){

    do {
      nvtxs  = graph->nvtxs;      
      map    = graph->cmap;
      cnvtxs = 0;           

      U_serializeIntArray(map,nvtxs,fp);
      graph = graph->coarser;
    }while(graph != NULL && graph->coarser != NULL);
  }

  U_serializeIntArray(cgraph->where,cgraph->nvtxs,fp);  
  
  fclose(fp);
  
}
/**
 * \brief Serialize the partition
 *
 * the partition is appended into [F_ouput].where.[id]
 *
 * \param graph partition to serialize
 * \param F_output IoCtrl::F_output
 * \param outid IoCtrl::out [IO_ID]
 * \param outresult IoCtrl::out [IO_RESULT]
 * \param outser IoCtrl::out [IO_SER]
 */
void Metis_SerializeWhere(GraphType *g,char* F_output,int outid,int outresult,int outser){
  FILE* fp;
  char out[256];
  if(outresult != IO_OUT_NO && (outser == IO_SER_DO || outser == IO_SER_USE)){
    if(outresult == IO_OUT_LEVEL)
      sprintf(out,"%s.where.%d",F_output,outid);
    else
      sprintf(out,"%s.where",F_output);
    NOTICE(((fp = fopen(out,"a+b")) == NULL,"MD_PagridWriteWhere error code 1\n"));
    U_serializeIntArray(g->where,g->nvtxs,fp);
    fclose(fp);
  }
}

/**
 * \brief unserialize a sequence of graph (matchings)
 *
 * \param F_ser file to unserialize
 * \param graph initial graph (if it was already readden) can be NULL
 * \return F_graph path to the graph file
 * \return F_grid path to the grid file
 * \return level number of level in the coarsening phase
 * \return CType coarsening scheme IoCtrl::CType
 * \return WType IoCtrl::WType
 * \return RType IoCtrl::RType
 * \return CSeed IoCtrl::CSeed
 * \return ISeed IoCtrl::ISeed
 * \return RSeed IoCtrl::RSeed
 *
 * \return Coarsest graph
 */
GraphType *Metis_Unserialize(char *F_ser,GraphType *graph,
			      char **F_graph,char **F_grid,
			      int* level, 
			      int* CType, int* WType, int *RType,
			      int* CSeed, int *ISeed, int *RSeed){
  int i;
  FILE *fp, *fpgraph;
  GraphType *fgraph;

  fp = PLIB_UnserializeInit(&i,F_ser, F_graph, F_grid,
			    level, CType, WType, RType,
			    CSeed, ISeed, RSeed);

  WARNING((i != TOOL_METIS,"Metis_Unserialize error code 1 : file %s was serialized by %s\n",F_ser,(i == TOOL_PAGRID?"PaGRID":"unknown tool")),NULL);
 
  Metis_Randini(*CSeed);

  if(graph == NULL){
    WARNING(((fpgraph = fopen(*F_graph,"r")) == NULL,"MD_PagridReadSerial error code 2 : cannot open %s\n",*F_graph),NULL);
    graph = __Read(fp);
    fclose(fpgraph);
  }

  fgraph = graph;

  for(i=0;i<*level;i++){
    U_unserializeIntArray(graph->nvtxs,fp,graph->cmap,LEND);
    graph       = __Buildcoarse(graph);
  }

  graph->where = U_unserializeIntArray(graph->nvtxs,fp,LEND);

  fclose(fp);

  return graph;
}

GraphType* __Buildcoarse(GraphType* graph){
    int i, idx, j, v, u, k, maxdegree=0;
  
  int fnvtxs     = graph->nvtxs;
  int fnedges    = graph->nedges;
  int* vwgt      = graph->vwgt;
  int* adjidx    = graph->xadj;
  int* adjncy    = graph->adjncy;
  int* adjwgt    = graph->adjwgt;
  int* adjwgtsum = graph->adjwgtsum;
  int* map       = graph->cmap;
  int* match     = U_malloc(fnvtxs,int,"Pagrid_Buildcoarse error code 1\n");
  int* mapidx    = U_mallocset(fnvtxs,int,-1,"Pagrid_Buildcoarse error code 1\n");

  GraphType *child;
  int *mapadj;
  int *record;

  int 
    nedges,  ncon,
    cnedges, cnvtxs = 0;
  
  int 
    *cadjncy, *cadjwgt, *cadjwgtsum, *cvwgt, *cadjidx;  

  /* UNW */
  if(!vwgt){
    vwgt      = U_mallocset(fnvtxs,int,1,"Metis_TransGraph error code 1\n");
    adjwgt    = U_mallocset(fnedges,int,1,"Metis_TransGraph error code 2\n");
    adjwgtsum = U_malloc(fnvtxs,int,"Metis_TransGraph error code 3\n");
    for(i=0;i<fnvtxs;i++){
      adjwgtsum[i] = adjidx[i+1] - adjidx[i];
    }
  }

  /* rebuild match */
  for(i=0;i<fnvtxs;i++){
    idx = map[i];

    if(idx+1>cnvtxs)
      cnvtxs = idx+1;

    if(mapidx[idx] == -1)
      mapidx[idx] = i;
    else{
      match[i] = mapidx[idx];
      match[mapidx[idx]] = i;
      mapidx[idx] = -2;
    }
  }

  for(i=0;i<fnvtxs;i++){
    idx = map[i];
    if(mapidx[idx] != -2){
      mapidx[idx] = -2;
      match[i] = i;
    }
  }

  /* rebuild guide */
  Metis_Perm(fnvtxs, mapidx);
  
  /* contract */
  child         = __Init();
  child->nvtxs  = cnvtxs;
  ncon          = child->ncon  = graph->ncon;
  child->nedges = fnedges;

  __Allocate(child);

  cadjidx    = child->xadj;
  cadjncy    = child->adjncy;
  cvwgt      = child->vwgt;
  cadjwgt    = child->adjwgt; 
  cadjwgtsum = child->adjwgtsum;
  
  cnedges = cadjidx[0] = 0;

  /*Initialize the vertices list of the coarser graph (-1 => not contracted)*/
  for (i =0; i<cnvtxs; i++)
    child->vwgt[i] = -1;
  
  record  = U_mallocset(cnvtxs,int,-1,"Contract error code 2\n");
  mapadj  = U_malloc(fnedges,int,"Contract error code 3\n");

  /* map the adjacency matrix of the finer graph */
  for (i = 0; i < fnedges; i++)
    mapadj[i] = map[adjncy[i]];

  /*****************************/
  /* Contraction               */
  /*****************************/
  cnvtxs = 0;
  for (i = 0; i < fnvtxs; i++) {
    v = mapidx[i];
    
    /*if the label to which vertex v is mapped is equal to cnvtxs, 
      collapse it and the vertex to which it matches*/
    if (map[v] != cnvtxs)
      continue;
      
    if (ncon == 1)
      cvwgt[cnvtxs] = vwgt[v];

    /*****************************/
    /* map v in coarser graph    */
    /*****************************/
    cadjwgtsum[cnvtxs] = adjwgtsum[v];	/*Initilize adjwgsum of 
					  vertex cnvtxs in coarser 
					  graph to the adjwgtsum of 
					  vertex v in the father graph*/
    
    nedges = 0;	/*Initialize number of edges of map[v] to be zero*/
    
    /* add all edges of v in the finer graph */
    for (j = adjidx[v]; j < adjidx[v+1]; j++) {
      k = mapadj[j];	
      /* if edges is not already existing add it */
      if (record[k] == -1) {
	cadjncy[nedges] = k;
	cadjwgt[nedges] = adjwgt[j];
	record[k] = nedges++;
      /* else only update its weight */
      } else {
	cadjwgt[record[k]] += adjwgt[j];
      }
    }
    
    /**********************************************/
    /* map u (matched vertex) in coarser graph    */
    /**********************************************/
    u = match[v];	
    if (v != u) {			
      if (ncon == 1)
	cvwgt[cnvtxs] += vwgt[u];	
      
      cadjwgtsum[cnvtxs] += adjwgtsum[u];	
      
      /* add all edges of u in the finer graph */
      for (j = adjidx[u]; j < adjidx[u+1]; j++) {
	k = mapadj[j];
	/* if edges is not already existing add it */
	if (record[k] == -1) {
	  cadjncy[nedges] = k;
	  cadjwgt[nedges] = adjwgt[j];
	  record[k] = nedges++;
	/* else only update its weight */
	} else {
	  cadjwgt[record[k]] += adjwgt[j];
	}
      }
      
      /* Remove the contracted edge if existing */
      if ((j = record[cnvtxs]) != -1) {
	cadjwgtsum[cnvtxs] -= cadjwgt[j];	/* subtract its weight */
	cadjncy[j] = cadjncy[--nedges];		/* replace with last edge */
	cadjwgt[j] = cadjwgt[nedges];		
	record[cnvtxs] = -1;
      }
    }
    /*********************/
    /* clean and update  */
    /*********************/
    for (j=0; j<nedges; j++)
      record[cadjncy[j]] = -1;	
    
    cnedges += nedges;			/*add up the number of edges 
					  found in a iteration*/	
    cadjidx[++cnvtxs] = cnedges;	/*add up the edges found in a 
					  iteration to adjidx of the 
					  coarser graph*/
    
    cadjncy += nedges;		/*move the pointer preparing 
				  for receving new adjacent edges*/
    cadjwgt += nedges;		/*move the pointer preparing 
				  for receving new adjacent edges*/
    if(nedges>maxdegree)
      maxdegree = nedges;
  }
  /************/
  /* finalize */
  /************/
  child->nedges = cnedges;

  /*attach child to graph*/
  graph->coarser = child;
  child->finer   = graph;
  
  /*free tmp arrays*/
  U_free("Contract error code 3\n",&match,&mapidx,&mapadj,&record,LEND);

  return child;
}

#endif
/** 
 * \brief return a string describing some of the options used
 *
 * The returned string must be freed
 * \param CType matching scheme
 * \param WType edge weight type during the coarsening phase
 * \param RType uncoarsening method
 * \param CSeed Coarsening seed
 * \param ISeed Initial partitioning Seed
 * \param RSeed Uncoarsening seed
 * \return string describing the tool used and some of its options
 */
char *Metis_GetScheme(int CType,int WType,int RType,int CSeed, int ISeed, int RSeed){
  char method[256];
  char *scheme;

  switch(CType){
  case IO_MATCH_LAM:
    strcpy(method , "Metis-LAM");
    break;
  case IO_MATCH_MWM:
    strcpy(method , "Metis-MWM");
    break;
  case IO_MATCH_RM:
    strcpy(method , "Metis-RandM");
    break;
  case IO_MATCH_HEM:
    strcpy(method , "Metis-HEM");
    break;
  case IO_MATCH_SHEM:
    strcpy(method , "Metis-MHEM");
    break;
  case IO_MATCH_SHEMKWAY:
    strcpy(method , "Metis-SHEMKWAY");
    break;
  case IO_MATCH_SHEBM_ONENORM:
    strcpy(method , "Metis-SHEBM_ONENORM");
    break;
  case IO_MATCH_SHEBM_INFNORM:
    strcpy(method , "Metis-SHEBM_INFNORM");
    break;
  case IO_MATCH_SBHEM_ONENORM:
    strcpy(method , "Metis-SBHEM_ONENORM");
    break;
  case IO_MATCH_SBHEM_INFNORM:
    strcpy(method , "Metis-SBHEM_INFNORM");
    break;
  default:
    strcpy(method , "Metis-unknown");
  }
  
  switch(WType){
  case IO_WEIGHT_NORM:
    strcat(method , "-NORM");
    break;
  case IO_WEIGHT_UPGMA:
    strcat(method , "-UPGMA");
    break;
  default:
    strcat(method , "-unknown");
  }

  if(CSeed != -1 && CSeed != IO_SEED)
    strcat(method,"-RM");
  if(ISeed != -1 && ISeed != IO_SEED)
    strcat(method,"-RI");
  if(RSeed != -1 && RSeed != IO_SEED)
    strcat(method,"-RR");
  
  scheme = U_malloc(strlen(method)+1,char,"MD_MetisGetScheme error code 1\n");
  strcpy(scheme,method);
  return scheme;
}
/**
 * \brief permute an array
 * \params n size of the array
 * \params arr arry to be permuted
 */
void Metis_Perm(int n,int* p){
  int i, u, v;
  idxtype tmp;
  
  for (i=0; i<n; i++)
    p[i] = i;

  if (n <= 4)
    return;

  for (i=0; i<n; i+=16) {
    u = RandomInRangeFast(n-4);
    v = RandomInRangeFast(n-4);
    SWAP(p[v], p[u], tmp);
    SWAP(p[v+1], p[u+1], tmp);
    SWAP(p[v+2], p[u+2], tmp);
    SWAP(p[v+3], p[u+3], tmp);
  }
}
/**
 * \brief initialize the random number generator
 * \params seed random number seed
 */
void Metis_Randini(int seed){
  if (seed == -1) {
    srand(4321);  
  }
  else {
    srand(seed);  
  }
}
